{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f1 的 LaTeX 代码:\n",
      "- 240.0 \\nu x^{4} y + 120.0 \\nu x^{4} + 480.0 \\nu x^{3} y - 240.0 \\nu x^{3} - 480.0 \\nu x^{2} y^{3} + 720.0 \\nu x^{2} y^{2} - 480 \\nu x^{2} y + 120.0 \\nu x^{2} + 480.0 \\nu x y^{3} - 720.0 \\nu x y^{2} + 240.0 \\nu x y - 80.0 \\nu y^{3} + 120.0 \\nu y^{2} - 40.0 \\nu y + 40 y - 20\n",
      "\n",
      "f2 的 LaTeX 代码:\n",
      "480.0 \\nu x^{3} y^{2} - 480.0 \\nu x^{3} y + 80.0 \\nu x^{3} - 720.0 \\nu x^{2} y^{2} + 720.0 \\nu x^{2} y - 120.0 \\nu x^{2} + 240.0 \\nu x y^{4} - 480.0 \\nu x y^{3} + 480 \\nu x y^{2} - 240.0 \\nu x y + 40.0 \\nu x - 120.0 \\nu y^{4} + 240.0 \\nu y^{3} - 120.0 \\nu y^{2} + 40 x - 20\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "import numpy as np\n",
    "\n",
    "# 定义符号变量\n",
    "x, y, nu = sp.symbols('x y nu')\n",
    "\n",
    "# 定义速度场 u 和压力场 p\n",
    "u1 = 20 * x**2 * (x - 1)**2 * y * (y - 1) * (2*y - 1)\n",
    "u2 = -20 * x * (x - 1) * (2*x - 1) * y**2 * (y - 1)**2\n",
    "p = 10 * (2*x - 1) * (2*y - 1)\n",
    "\n",
    "# 定义 D(u)\n",
    "\n",
    "\n",
    "def D(u):\n",
    "    du1_dx = sp.diff(u[0], x)\n",
    "    du1_dy = sp.diff(u[0], y)\n",
    "    du2_dx = sp.diff(u[1], x)\n",
    "    du2_dy = sp.diff(u[1], y)\n",
    "    return [[du1_dx, 0.5*(du1_dy + du2_dx)], [0.5*(du1_dy + du2_dx), du2_dy]]\n",
    "\n",
    "# 定义 T(u, p)\n",
    "\n",
    "\n",
    "def T(u, p):\n",
    "    D_u = D(u)\n",
    "    return [[2*nu*D_u[0][0] - p, 2*nu*D_u[0][1]], [2*nu*D_u[1][0], 2*nu*D_u[1][1] - p]]\n",
    "\n",
    "\n",
    "# 定义速度场 u 为向量\n",
    "u = [u1, u2]\n",
    "\n",
    "# 计算 T(u, p)\n",
    "T_u_p = T(u, p)\n",
    "\n",
    "# 计算 T(u, p) 的散度\n",
    "div_T_x = sp.diff(T_u_p[0][0], x) + sp.diff(T_u_p[0][1], y)\n",
    "div_T_y = sp.diff(T_u_p[1][0], x) + sp.diff(T_u_p[1][1], y)\n",
    "\n",
    "# 计算体积力 f\n",
    "f1 = -div_T_x\n",
    "f2 = -div_T_y\n",
    "\n",
    "# 整理表达式\n",
    "f1_expanded = sp.expand(f1)\n",
    "f1_simplified = sp.simplify(f1_expanded)\n",
    "f2_expanded = sp.expand(f2)\n",
    "f2_simplified = sp.simplify(f2_expanded)\n",
    "\n",
    "\n",
    "# 输出 f1 和 f2 的简化形式的 LaTeX 代码\n",
    "print(\"f1 的 LaTeX 代码:\")\n",
    "print(sp.latex(f1_simplified))\n",
    "print(\"\\nf2 的 LaTeX 代码:\")\n",
    "print(sp.latex(f2_simplified))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dolfinx",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
